correction of the modulated binning

In [1]:

    

import PyPlot
import DSP


    

In [2]:

    

include("../juwvid.jl")


    

    
Out[2]:





juwvid

In [3]:

    

## Generating the modulated time bins (td) from the unifrom time bins (t)
nsample=4024;
dt=30.0/60/24
t=collect(1:nsample)*dt;

p=1.e-3 #modulation factor
wbin=1/5.0
td=t.+p*sin.(wbin*t);


    

In [4]:

    

PyPlot.plot(t,t-td,".")
PyPlot.ylabel("time - modulated time")
PyPlot.xlabel("time")


    

    













    
Out[4]:





PyObject Text(0.5,24,'time')

In [5]:

    

wp=2*pi
y=cos.(wp*td); # modulated signal
z=DSP.Util.hilbert(y);

yt=cos.(wp*t); # non-modulated signal for comparison
zt=DSP.Util.hilbert(yt);


    

In [6]:

    

fig = PyPlot.figure(figsize=(10,3))
PyPlot.plot(td,y,".",alpha=0.5)
PyPlot.plot(td,yt,".",alpha=0.5)
PyPlot.xlim(0,30)
PyPlot.ylim(-1.2,1.2)


    

    













    
Out[6]:





(-1.2, 1.2)

In [7]:

    

tfrpf=cohenclass.tfrpwv(z,NaN,NaN,NaN,NaN,NaN,0);


    

    




Single pseudo Wigner Ville
Use fft.

In [8]:

    

fig=PyPlot.figure()
ax = fig[:add_subplot](1,2,1)
a=juwplot.tfrshow(abs.(tfrpf),td[2]-td[1],td[1],td[end],NaN,NaN,0.7,"CMRmap")
PyPlot.xlabel("time [day]")
PyPlot.ylabel("frequency [1/day]")
PyPlot.title("Pseudo WV")
ax = fig[:add_subplot](1,2,2)
a=juwplot.tfrshow(abs.(tfrpf),td[2]-td[1],td[1],td[end],NaN,NaN,0.7*560,"CMRmap")
PyPlot.xlabel("time [day]")
PyPlot.title("Zoom-up")
PyPlot.ylim(0.98,1.02)


    

    













    
Out[8]:





(0.98, 1.02)

In [9]:

    

nft=512


    

    
Out[9]:




512

In [10]:

    

js,je=juwutils.frequency_to_index([0.999,1.001], td[2]-td[1], nsample, nft)


    

    
Out[10]:





2-element Array{Float64,1}:
 167.532
 167.868

In [11]:

    

#check
juwutils.index_to_frequency([js,je],NaN,td[2]-td[1], nsample, nft, je,js)


    

    
Out[11]:





2-element Array{Float64,1}:
 0.999
 1.001

In [12]:

    

nft=1024
fin=collect(linspace(js,je,nft));
tfrpfn=cohenclass.tfrpwv(z,NaN,NaN,fin,NaN,NaN,0,"nufft",16);
tfrpfnt=cohenclass.tfrpwv(zt,NaN,NaN,fin,NaN,NaN,0,"nufft",16);


    

    




Single pseudo Wigner Ville
Use nufft.
Single pseudo Wigner Ville
Use nufft.

In [13]:

    

indfn=extif.maxif(abs.(tfrpfn));
fn=juwutils.index_to_frequency(indfn,fin,t[2]-t[1],nsample,nft);
indfnt=extif.maxif(abs.(tfrpfnt));
fnt=juwutils.index_to_frequency(indfnt,fin,t[2]-t[1],nsample,nft);


    

In [14]:

    

# computing teh correction factor
dtddt=Float64[]
for i=1:length(t)-1
    de=(td[i+1]-td[i])/dt
    #println(de)
    append!(dtddt,[de])
end

PyPlot.plot(t[1:end-1],dtddt)
PyPlot.title("Correction Factor")


    

    













    
Out[14]:





PyObject Text(0.5,1,'Correction Factor')

In [15]:

    

fig=PyPlot.figure(figsize=(10,3))
ax = fig[:add_subplot](1,2,1)
PyPlot.xlabel("time")
PyPlot.ylabel("instantaneous frequency")
PyPlot.ylim(0.999,1.001)
PyPlot.plot(td,fn, color="C0",lw=1)
PyPlot.plot(t[1:end-1],dtddt, color="C3",lw=2,ls="dashed")
PyPlot.plot(t[1:end-1],fn[1:end-1]./dtddt, color="C2",lw=1)
PyPlot.legend(["modulated signal","expected modulation","after correction"])

ax = fig[:add_subplot](1,2,2)
PyPlot.plot(t,fnt, color="C1",ls="dashed",lw=1)
PyPlot.xlabel("time")
#PyPlot.ylabel("instantaneous frequency")
PyPlot.ylim(0.999,1.001)
PyPlot.legend(["original signal"])


    

    













    
Out[15]:





PyObject <matplotlib.legend.Legend object at 0x7f0ac934c978>

In [ ]: